dir = pwd();
dir in LOAD_PATH ? nothing : push!(LOAD_PATH, dir);
println(LOAD_PATH[end])
using FEM_module, Dqcoeff, OffsetArrays, PyCall, ProgressMeter
plt = pyimport("matplotlib.pyplot")
np = pyimport("numpy")
include("parameters.jl")
include("creep_and_shrinkage_functions.jl")
include("shape_functions.jl")
include("ke_loade.jl")

nels = 12;  		#(*total element number*)
ends = subnodes(Int64(nels/2), 2.0, 0.0, span);
nodof = 10; 		#(*dof number of each node*)
nod = ndp;    		#(*node number of each element*)
ndof = nod*nodof; 	#(*dof number of each element*)

gnum = zeros(Int64, nels, nod)
for iel = 1:nels
	gnum[iel,:] = range( 1+(iel-1)*(ndp-1), stop = ndp+(iel-1)*(ndp-1) )
end

elcd = []
for iel = 1:nels
	v = 0.5*(ends[iel + 1] - ends[iel])*expanded_chebyshev(ndp) .+ 0.5*(ends[iel + 1] + ends[iel])
	if iel == 1
		global elcd = [elcd; v]
	else
		global elcd = [elcd; v[2:end]]
	end
end

nn = length(elcd) # total node number
gcoord = Array{Float64,2}(reshape(elcd,nn,1))

nf = ones(Int64, nodof, nn)
nr = 1; # number of dof which is suppressed
fix = [1 0 0 0 0 0 0 0 0 0 0;
      nn 0 0 0 0 0 0 0 0 0 0;]
for i = 1:nr
	nf[:, fix[i, 1]] = fix[i, 2 : end]
end
formnf!(nf)
neq = maximum(nf)
g_g = gg(nf, ndof, gnum)
kdiag = zeros(Int64, neq)
for iel = 1:nels
	g = g_g[: , iel]
	fkdiag!(kdiag, g)
end
for i = 2:neq
	kdiag[i] = kdiag[i] + kdiag[i - 1]
end

sols = zeros(length(time0), neq)
StressHistory = zeros(length(time0), nels, 3*ndp)
StressHistorytao = zeros(length(time0), nels, 3*ndp)

# global stiffness matrix assembly
kv0 = zeros(kdiag[neq])
for iel = 1:nels
	num = gnum[iel,:]
	coord = gcoord[num,:]
	(xi, xj) = coord[[1,ndp],1]
	km0 = ksteel(xi, xj) +  kcs(xi,xj)
	g = g_g[: , iel]
	fsparv!(kv0, km0, g, kdiag)
end
loads0 = assemble_loads0(gnum, gcoord, nf, nels, neq)

@showprogress 1.0 "computing... " for local_nt = 1:length(time0)
	global nt = local_nt
	global Ec = en(nt,time0)*1e6
	global Gc = beta_E(time0[nt])*E28/(2*(1 + μc))
	kv1 = zeros(kdiag[neq])
	for iel = 1:nels
		num = gnum[iel,:]
		coord = gcoord[num,:]
		(xi, xj) = coord[[1,ndp],1]
		km1 = kconcrete(xi, xj) + krebar(xi, xj)
		g = g_g[: , iel]
		fsparv!(kv1, km1, g, kdiag)
	end
	loads1 = assemble_loads1(gnum, gcoord, nf, nels, neq)

	kv = kv0 + kv1; loads = loads0 + loads1;
	sparin_gauss!(kv, kdiag)
	spabac_gauss!(kv, loads, kdiag)
	sols[nt,:] = loads
	record(loads, nt)
end

# ==================== post_processing =====================
include("post_processing.jl")
